library(tidyverse)
library(gifski)
library(ggraph)
library(here)
library(igraph)
source(here("modelFunction_rewiring.R"))
First, run the model.
# Define parameters
N = 100
edge.prob <- 0.04
burn.in = 20
burn.out = 5
pm = 0.3
ps = 0.1
pa = 0.2
add00 = c(0.5, 10)
lose01 = 0.1
add10 = 0.05
lose11 = c(0.5, 0.5)
histMultiplier = 1.2
doRemoval = TRUE
modelGraphs <- runModel(N = N, # Nodes in the network
edge.prob = edge.prob,
burn.in = burn.in,
burn.out = burn.out,
pm = pm,
ps = ps,
pa = pa,
add00 = add00,
lose01 = lose01,
add10 = add10,
lose11 = lose11,
histMultiplier = histMultiplier,
doRemoval = doRemoval) %>%
lapply(., function(x){
igraph::graph_from_adjacency_matrix(x, mode = "undirected", add.colnames = "label")
})
Compute degree and betweenness for each of the model networks.
attrDF <- lapply(modelGraphs, function(x){
x %>%
set_vertex_attr(., name = "degree", value = degree(x)) %>%
set_vertex_attr(., name = "betweenness", value = betweenness(x)) %>%
tidygraph::as_tbl_graph() %>%
tidygraph::activate(nodes) %>%
data.frame()
}) %>%
data.table::rbindlist(idcol = "slice") %>%
as.data.frame() %>%
mutate(beforeAfter = ifelse(slice <= burn.in, "before", "after"))
# make a long-format version to use for facetted plots
attrDFLong <- attrDF %>%
pivot_longer(cols = c(degree, betweenness), names_to = "measure", values_to = "value")
View degree and betweenness over time
# Degree and betweenness over time
attrDFLong %>%
ggplot(aes(x = slice, y = value, col = label, group = label))+
facet_wrap(~measure, scales = "free")+
geom_line(alpha = 0.7, size = 0.5)+
theme_minimal()+
geom_vline(xintercept = burn.in+1, col = "red", size = 1)+
theme(legend.position = "none")
# Let's just look at a few labels, because this plot is too hard to read
indivs <- sample(unique(attrDFLong$label), 5)
attrDFLong %>%
filter(label %in% indivs) %>%
ggplot(aes(x = slice, y = value, col = label, group = label))+
facet_wrap(~measure, scales = "free")+
geom_line()+
theme_minimal()+
geom_vline(xintercept = burn.in+1, col = "red", size = 1)
# This is still a little hard to see... let's look at just the "before" rows.
attrDFLong %>%
filter(label %in% indivs, beforeAfter == "before") %>%
ggplot(aes(x = slice, y = value, col = label, group = label))+
facet_wrap(~measure, scales = "free")+
geom_line()+
theme_minimal()
# Okay, I'm seeing some differentiation, but no clear patterns. Because of how the model is set up, the rank order tends to be maintained for a few consecutive time steps, but not over the whole span of the baseline model. This makes sense, since labels don't have any a priori reason to have higher or lower degrees than other labels.
Now, going to compute three network-level measures before and after the removal of an label.
How does edge density behave over time?
densities <- lapply(modelGraphs, function(x){
edge_density(x)
}) %>%
unlist() %>%
as.data.frame() %>%
setNames(., "density") %>%
mutate(slice = 1:length(modelGraphs))
densities %>%
ggplot(aes(x = slice, y = density))+
geom_line()+
theme_minimal()+
geom_vline(aes(xintercept = burn.in+1), size = 1, col = "red")
First delta: immediately before removal to after removal/pre-rewiring
delta1_density <- igraph::edge_density(modelGraphs[[burn.in+1]])-igraph::edge_density(modelGraphs[[burn.in]])
Second delta: after removal/pre-rewiring to post-rewiring
delta2_density <- igraph::edge_density(modelGraphs[[burn.in+2]])-igraph::edge_density(modelGraphs[[burn.in+1]])
How does connectivity behave over time?
connectivities <- lapply(modelGraphs, function(x){
mean_distance(x)
}) %>%
unlist() %>%
as.data.frame() %>%
setNames(., "connectivity") %>%
mutate(slice = 1:length(modelGraphs))
connectivities %>%
ggplot(aes(x = slice, y = connectivity))+
geom_line()+
theme_minimal()+
geom_vline(aes(xintercept = burn.in+1), size = 1, col = "red")
## Warning: Removed 1 row(s) containing missing values (geom_path).
First delta: immediately before removal to after removal/pre-rewiring
delta1_meanDistance <- igraph::mean_distance(modelGraphs[[burn.in+1]])-igraph::mean_distance(modelGraphs[[burn.in]])
Second delta: after removal/pre-rewiring to post-rewiring
delta2_meanDistance <- igraph::mean_distance(modelGraphs[[burn.in+2]])-igraph::mean_distance(modelGraphs[[burn.in+1]])
How does modularity behave over time?
clustered <- lapply(modelGraphs, function(x){
x %>%
cluster_fast_greedy()
})
modularities <- data.frame(slice = 1:length(clustered),
modularity = unlist(lapply(clustered, modularity)),
nClusters = unlist(lapply(clustered, length)),
nIndivs = unlist(lapply(clustered, function(x){length(membership(x))})))
modularities %>%
ggplot(aes(x = slice, y = modularity))+
geom_line()+
theme_minimal()+
geom_vline(aes(xintercept = burn.in+1), size = 1, col = "red")
modularities %>%
ggplot(aes(x = slice, y = nClusters))+
geom_line()+
theme_minimal()+
geom_vline(aes(xintercept = burn.in+1), size = 1, col = "red")
modularities %>%
ggplot(aes(x = slice, y = nClusters/modularity))+
geom_line()+
theme_minimal()+
geom_vline(aes(xintercept = burn.in+1), size = 1, col = "red")
# there doesn't seem to be a pattern in the ratio of nClusters to modularity, but I will investigate on a larger scale.
First delta: immediately before removal to after removal/pre-rewiring
delta1_modularity <- modularity(cluster_fast_greedy(modelGraphs[[burn.in+1]]))-modularity(cluster_fast_greedy(modelGraphs[[burn.in]]))
Second delta: after removal/pre-rewiring to post-rewiring
delta2_modularity <- modularity(cluster_fast_greedy(modelGraphs[[burn.in+2]]))-modularity(cluster_fast_greedy(modelGraphs[[burn.in+1]]))
I am not at all convinced that these momentary deltas are what I should be using. Seems like we need to at least consider two time steps out, since the model relies on this. Looking at the graphs above shows me how different the longer-term trends are from the short-term fluctuations.
Let’s run the model 100 times and see if trends emerge in how to measure.
nRuns <- 100
modelRuns <- vector(mode = "list", length = nRuns)
for(i in 1:nRuns){
modelRuns[[i]] <- runModel(N = N, # Nodes in the network
edge.prob = edge.prob,
burn.in = burn.in,
burn.out = burn.out,
pm = pm,
ps = ps,
pa = pa,
add00 = add00,
lose01 = lose01,
add10 = add10,
lose11 = lose11,
histMultiplier = histMultiplier,
doRemoval = doRemoval) %>%
lapply(., function(x){
igraph::graph_from_adjacency_matrix(x, mode = "undirected")
})
}
densList <- lapply(modelRuns, function(x){
lapply(x, function(y){
edge_density(y)
}) %>%
unlist() %>%
as.data.frame() %>%
setNames(., "density") %>%
mutate(slice = 1:nrow(.))
})
densDF <- data.table::rbindlist(densList, idcol = "run")
connList <- lapply(modelRuns, function(x){
lapply(x, function(y){
mean_distance(y)
}) %>%
unlist() %>%
as.data.frame() %>%
setNames(., "connectivity") %>%
mutate(slice = 1:nrow(.))
})
connDF <- data.table::rbindlist(connList, idcol = "run")
moduList <- lapply(modelRuns, function(x){
clustered <- lapply(x, cluster_fast_greedy)
df <- data.frame(slice = 1:length(clustered),
modularity = unlist(lapply(clustered, modularity)),
nClusters = unlist(lapply(clustered, length)),
nIndivs = unlist(lapply(clustered, function(i){length(membership(i))})))
return(df)
})
moduDF <- data.table::rbindlist(moduList, idcol = "run")
Make some plots and try to detect any changepoint trends
densDF %>%
ggplot(aes(x = slice, y = density))+
geom_line(aes(group = run, col = run), size = 0.1)+
theme_minimal()+
geom_vline(aes(xintercept = burn.in+1), size = 0.5, col = "red")+
scale_color_viridis_c()
connDF %>%
ggplot(aes(x = slice, y = connectivity))+
geom_line(aes(group = run, col = run), size = 0.1)+
theme_minimal()+
geom_vline(aes(xintercept = burn.in+1), size = 0.5, col = "red")+
scale_color_viridis_c()
## Warning: Removed 100 row(s) containing missing values (geom_path).
moduDF %>%
ggplot(aes(x = slice, y = modularity))+
geom_line(aes(group = run, col = run), size = 0.1)+
theme_minimal()+
geom_vline(aes(xintercept = burn.in+1), size = 0.5, col = "red")+
scale_color_viridis_c()
moduDF %>%
filter(slice > 3) %>%
ggplot(aes(x = slice, y = nClusters))+
geom_line(aes(group = run, col = run), size = 0.1)+
theme_minimal()+
geom_vline(aes(xintercept = burn.in+1), size = 0.5, col = "red")+
scale_color_viridis_c()
moduDF %>%
ggplot(aes(x = slice, y = modularity/nClusters))+
geom_line(aes(group = run, col = run), size = 0.1)+
theme_minimal()+
geom_vline(aes(xintercept = burn.in+1), size = 0.5, col = "red")+
scale_color_viridis_c()
# Absolutely no pattern at all in the ratio of clusters to modularity.